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ABSTRACT 

The  focus  of  this  paper  is  on  the  development  of  methodology  for  short-term  (1-3  days)  oceanic  biolu¬ 
minescence  (BL)  predictions  and  the  optimization  of  spatial  and  temporal  bioluminescence  sampling  strat¬ 
egies.  The  approach  is  based  on  predictions  of  bioluminescence  with  an  advection-diffusion-reaction 
(tracer)  model  with  velocities  and  diffusivities  from  a  circulation  model.  In  previous  research,  it  was  shown 
that  short-term  changes  in  some  of  the  salient  features  in  coastal  bioluminescence  can  be  explained  and 
predicted  by  using  this  approach.  At  the  same  time,  it  was  demonstrated  that  optimization  of  biolumines- 
ccncc  sampling  prior  to  the  forecast  is  critical  for  successful  short-term  BL  predictions  with  the  tracer 
model.  In  the  present  paper,  the  adjoint  to  the  tracer  model  is  used  to  study  the  sensitivity  of  the  modeled 
bioluminescence  distributions  to  the  sampling  strategies  for  BL.  The  locations  and  times  of  bioluminescence 
sampling  prior  to  the  forecast  are  determined  by  using  the  adjoint-based  sensitivity  maps.  The  approach  is 
tested  with  bioluminescence  observations  collected  during  August  2000  and  2003  in  the  Monterey  Bay, 
California,  area.  During  August  2000,  BL  surveys  were  collected  during  a  strong  wind  relaxation  event, 
while  in  August  2003.  BL  surveys  were  conducted  during  an  extended  (longer  than  a  week)  upwelling- 
favorable  event.  The  numerical  bioluminescence  predictability  experiments  demonstrated  a  close  agree¬ 
ment  between  observed  and  model-predicted  short-term  spatial  and  temporal  changes  of  the  coastal  bio¬ 
luminescence. 


1.  Introduction 

In  previous  research  (Shulman  et  al.  2003),  it  was 
demonstrated  that  some  of  the  salient  short-term 
changes  in  bioluminescence  (BL)  intensity  can  be 
explained  by  hydrodynamic  transport  processes.  Bio- 


*  Nava  I  Research  Laboratory  Contribution  Number  JA/7330- 
04-10,  AOSN  Contribution  Number  2005.101,  and  Woods  Hole 
Oceanographic  Institution  Contribution  Number  11253. 

Corresponding  author  address:  Igor  Shulman.  Naval  Research 
Laboratory.  Code  7331,  Bldg.  1009.  Stennis  Space  Center,  MS 
39529-5004. 

E-mail:  igor.shulmantohirlssc.na vy.mil 


luminescence  short-term  predictability  experiments 
were  conducted  by  assimilating  BL  observations  in¬ 
to  an  advection-diffusion-reaction  (tracer)  model 
with  velocities  and  diffusivities  from  a  circulation 
model. 

The  methodology  was  tested  with  BL  observations 
acquired  during  the  August  2000  experiment  conducted 
jointly  by  the  Autonomous  Ocean  Sampling  Network 
(AOSN  I),  Monterey  Bay  Aquarium  Research  Institute 
(MBARl)  Upper-Water-Column  Science  Experiment 
(MUSE),  and  National  Oceanic  Partnership  Program 
(NOPP)  Innovative  Coastal-Ocean  Observing  Network 
(ICON)  projects.  The  results  of  numerical  experiments 
designed  to  estimate  the  limits  of  bioluminescence  pre¬ 
dictions  by  the  tracer  model  demonstrated  the  strong 
utility  of  the  proposed  methodology  in  prediction  and 
interpretation  of  observed  short-term  changes  (1-3 
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days)  in  bioluminescence  intensity  (Shulman  et  al. 
2003). 

At  the  same  time,  results  show  that  optimization  of 
locations  and  times  of  BL  observations  prior  to  the 
forecast  are  critical  for  successful  short-term  BL  pre¬ 
dictions.  It  was  demonstrated  that  during  a  strong  wind 
relaxation  event  in  Monterey  Bay,  California,  the  as¬ 
similation  of  only  one  section  of  observed  BL  inside  the 
bay  gave  a  good  reconstruction  of  the  location  and  the 
maximum  of  BL  observed  (not  assimilated)  outside  of 
the  bay.  Also,  it  was  shown  that  sampling  of  BL  outside 
the  bay  would  provide  little  information  for  short-term 
BL  predictability  inside  the  bay.  For  short-term  BL 
predictability,  the  sampling  of  BL  intensity  should  be 
done  in  particularly  flow-dependent,  “sensitive'1  re¬ 
gions. 

Here,  short-term  BL  forecast  sensitivities  to  the 
circulation  patterns,  as  well  as  to  BL  intensity  (1-3 
days  prior  to  the  forecast)  are  presented.  Sensitivity 
studies  were  conducted  by  using  the  adjoint  to  the 
advection-diffusion  reaction  model.  Sensitivity  maps 
illustrate  specific  areas  where  initial  conditions  are 
most  influential  on  the  forecast  and  where  sampling 
for  BL  intensity  is  critical.  Bioluminescence  sensi¬ 
tivity  maps  are  presented  for  two  extensive  BL  sur¬ 
veys  during  the  AOSN  I  (August  2000)  and  AOSN  II 
(August  2003)  experiments.  Sensitivity  studies  were 
conducted  for  two  major  circulation  regimes  in 
Monterey  Bay — upwelling  (August  2003)  and  relax¬ 
ation  (August  2000).  Analysis  of  the  sensitivity  maps 
for  August  2000  demonstrated  good  agreement  with 
BL  predictability  results  presented  in  Shulman  et  al. 
(2003).  New  BL  predictability  experiments  were  con¬ 
ducted  utilizing  sensitivity  maps  for  the  August  2003 
experiment. 

The  paper  has  the  following  structure.  In  section  2, 
circulation  models  of  the  Monterey  Bay  area  are  briefly 
presented.  Section  3  describes  the  observed  BL  data 
used  in  this  study.  Section  4  discusses  the  estimation  of 
sensitivity  maps  with  the  adjoint  and  the  optimization 
of  BL  sampling  with  sensitivity  maps.  Bioluminescence 
predictability  experiments  during  the  upwelling  event 
of  2003  are  presented  in  section  5,  with  conclusions  and 
future  plans  in  section  6. 

2.  Circulation  models 

Two  circulation  models  of  Monterey  Bay  are  used 
in  this  study  (Fig.  1).  The  ICON  model  has  an  orthogo¬ 
nal,  curvilinear  grid  with  horizontal  resolution  rang¬ 
ing  from  1  to  4  km.  The  model  has  30  vertical  sigma 


35. b  - 1 — 
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Fie.  1.  The  ICON  and  frslCON  model  domains  with  the 
locations  of  cross  sections  of  biolumincscence  observations. 


levels  (Shulman  el  al.  2002).  At  the  smaller  scale, 
a  fine-resolution  submodel  (frslCON)  has  been  set 
up  within  the  ICON  model  domain  (Fig.  1).  The  grid 
of  frslCON  has  a  variable  resolution  in  the  horizon¬ 
tal,  with  a  finer  resolution  (500-600  m)  around  the 
upwelling  front  in  the  northern  part  of  Monterey 
Bay  telescoping  to  a  coarser  resolution  (1.5  km)  in 
the  outer  portion  of  the  domain  (Shulman  et  al. 
2003).  The  ICON  and  frslCON  models  are  based 
on  the  three-dimensional,  sigma-coordinate  version 
of  the  Blumberg  and  Mellor  (1987)  hydrodynamic 
model. 

One  of  the  objectives  of  the  numerical  simula¬ 
tions  during  August  2000  was  the  modeling  of  the 
BL  intensity.  The  question  of  whether  short-term 
changes  in  some  features  of  BL  variability  can  be 
explained  by  hydrodynamic  transport  processes  was 
investigated.  For  this  reason,  the  focus  of  numerical 
simulations  was  on  accurate  model  predictions  in  the 
upper  100  m  of  the  water  column  (especially  predic¬ 
tions  of  velocity  fields).  To  achieve  this  objective,  the 
following  important  observations  and  forcings  were  uti¬ 
lized  during  the  August  2000  and  August  2003  simula¬ 
tions: 

(a)  Models  were  forced  with  high-resolution  (9  km 
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27  Aug  2000 


Surface 


Surface 


31  Aug  2000 


-122.4  -122.2  -122  -121.1 

Longitude 


50m 


50m 


Fig.  2.  ICON  model  surface  and  50-m-depth  currents  during  the  wind  relaxation  event  oi  27-31  Aug  2000. 


in  August  2000  and  3  km  in  August  2003)  wind 
stresses  and  heat  fluxes  from  the  Navy  Coupled 
Ocean  and  Atmospheric  Mesoscale  Prediction 
System  (COAMPS)  (Hodur  1997;  Kindle  et  al. 
2002). 

(b)  High-frequency  coastal  radar  (CODAR)  surface 
currents  were  assimilated  into  hydrodynamic  mod¬ 
els  based  on  the  scheme  described  in  Paduan  and 
Shulman  (2004). 

(c)  Open  boundary  conditions  for  the  ICON  model 
were  derived  from  the  larger-scale  Pacific  West 
Coast  (PWC)  model  predictions.  The  PWC  model 
domain  extends  seaward  to  135°W,  and  from  30°  to 
49°N.  The  PWC  model  was  forced  with  27-km 
COAMPS  wind  forcing,  and  the  multichannel  sea 
surface  temperature  (MCSST)  was  assimilated  into 
the  PWC  model. 

During  northwesterly,  upwelling  favorable  winds. 


the  hydrographic  conditions  in  and  around  Monterey 
Bay  are  mostly  determined  by  the  interaction  be¬ 
tween  upwelling  filaments  formed  at  headlands  to 
the  north  and  south  of  the  bay  and  anlicyclonic 
California  Current  meander  offshore  of  the  bay  [see 
Rosenfeld  et  al.  (1994)  for  more  details].  When 
upwelling-favorable  winds  weaken  (wind  relaxation), 
and  sometimes  become  poleward,  the  anlicyclonic  Cali¬ 
fornia  Current  meander  moves  onshore  and  then 
quickly  retreats  back  offshore  when  the  winds  reinten¬ 
sify. 

During  the  relaxation  event  of  August  2000  (Fig.  2), 
analysis  of  the  ICON  model  current  structure  at  vari¬ 
ous  depths  indicates  the  development  of  the  near¬ 
shore  northward  flow  extending  to  a  depth  of  50  m. 
This  northward  flow  is  connected  with  the  north¬ 
ward  flow  originating  at  the  southern  open  boundary 
of  the  ICON  domain  (a  larger-scale  phenomenon 
generated  by  the  coupling  to  the  PWC  model).  In 


1270 


JOURNAL  OF  ATMOSPHERIC  AND  OCEANIC  TECHNOLOGY 


VOLUMH  22 


Fig.  3.  ICON  model  surface  and  5()-m-depth  currents  during  the  upwelling  event  of  14-18  Aug  2003. 


the  bay,  this  northward  flow  develops  a  cyclonic 
circulation  that  is  confined  between  this  near¬ 
shore  northward  flow  and  the  southward  flow  off¬ 
shore. 

Analysis  of  the  ICON  model  current  structure  during 
the  upwelling  event  of  August  2003  (Fig.  3)  indicates  a 
strong  southward-flowing  offshore  jet.  This  southward 
jet  intensifies  and  flows  along  the  entrance  to  the  bay. 
At  the  same  time,  a  strong  cyclonic  eddy  is  present 
inside  the  bay. 


3.  Bioluminescence  observations 

Bioluminescence  data  in  Monterey  Bay  was  collected 
at  night  during  August  2003.  The  BL  was  measured 
using  two  custom-built  bathyphotometers,  each 
mounted  on  a  two  autonomous  underwater  vehicles 


(AUVs).  Both  AUVs  [a  Remote  Environmental  Mea¬ 
suring  Units  (REMUS)  (Moline  el  al.  2005)  and  a  Do¬ 
rado  (Wilcox  el  al.  2001)]  proceeded  through  the  water 
along  a  preprogrammed  path,  undulating  between  shal¬ 
low  (2  m)  and  deep  depth  boundaries  (40  m).  Within 
each  vehicle,  in  addition  to  the  core  instrument  pack¬ 
ages,  a  bioluminescence  bathyphotometer  (Elerren  et 
al.  2004)  pumps  water  into  a  0.5-L  sample  chamber  at  a 
rate  of  0.5  Ls-1.  Flow  rate,  temperature,  and  light  lev¬ 
els  (photon  flux,  assuming  isotropic  emission)  are  mea¬ 
sured  in  the  sample  chamber.  The  instrument  is  cali¬ 
brated  both  radiometrically  by  a  known  light  source 
and  biologically  by  insertion  of  a  known  concentration 
of  dinoflagellates.  As  with  all  bathyphotometers,  only  a 
fraction  of  the  total  luminescence  is  stimulated  or 
sampled,  but  because  organisms  emit  most  of  their  light 
in  the  first  part  of  their  first  flash,  we  capture  a  repeal- 
able  and  representative  fraction  of  the  bioluminescence 
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Fig.  4.  Observed  bioluminescence  distributions  along  cross  sections  in  Monterey  Bay  during  Aug  2003.  The  distance  scale  starts 

from  the  onshore  end  of  the  sections. 


present  in  the  environment.  This  component  is  the  bio¬ 
luminescence  that  we  are  modeling. 

Figure  4  shows  observed  BL  distributions  along  four 
sections.  Two  of  them  are  in  the  southern  portion  of  the 


bay  (see  map  on  Fig.  1):  the  section  labeled  H225  was 
taken  on  13  August,  while  FI226  was  taken  on  14  Au¬ 
gust.  The  other  two  sections  labeled  M14a  and  Ml4b  in 
the  northern  part  of  the  bay  were  taken  on  14  August. 
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Fig.  5.  Observed  bioluminescence  distributions  along  cross  section  M14a  during  15-18  Aug  2003.  The  distance  scale  starts  from  the 

onshore  end  of  the  section. 


The  observed  BL  distributions  (in  LO9  photons  per  sec¬ 
ond)  are  shown  as  a  function  of  depth  and  distance 
offshore. 

Another  set  of  BL  data  used  in  this  study  is  presented 


in  Fig.  5.  It  shows  the  spatial  and  temporal  change  in 
the  BL  intensity  over  4  days  (15-18  August  2003)  at 
section  M14a  (see  Fig.  1).  Figure  5  shows  offshore 
spreading  and  intensification  of  observed  BL  intensity 
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over  3-4  days  of  upwelling-favorable  conditions  in  the 
bay  (see  section  2  and  Fig.  3). 


4.  Sensitivity  studies  and  optimization  of 
bioluminescence  sampling 


In  Shulman  et  al.  (2003),  BL  predictability  experi¬ 
ments  were  conducted  with  the  use  of  the  advection- 
diffusion-reaction  equation 
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where  diffusivities  (A  and  K)  and  velocities  (u,  u,  vv)  are 
from  the  ICON  and  frslCON  circulation  models,  and 
S(x,  y,  z>  0  is  the  source-minus-sink  term  for  C.  In  this 
case  BL  is  modeled  as  concentration  C.  For  initializa¬ 
tion,  available  BL  observations  are  constantly  assimi¬ 
lated  into  the  model  (1)  by  using  the  source  term  S(x y} 
z>  t)  in  the  following  form  (Shulman  et  al.  2003): 

S(.t,  y,  Zs  t)  =  y(C  —  C°)8(t  -  r°),  (2) 

where  Cl)  are  BL  observations,  y  is  the  scalar  nudging 
coefficient  multiplying  (C  -  C°),  r  is  the  location  in  the 
model  domain  with  coordinates  (x,  y,  z),  r°  is  the  loca¬ 
tion  of  the  observed  BL  (C°)  with  coordinates  (x(\  y°9 
Z°),  and  8(t  -  r°)  is  a  Dirac  function  for  which  8  =  1 
when  t  =  t(1  and  8  =  0  for  all  other  cases. 

Velocities  and  diffusivities  in  (l)  are  taken  from  the 
initialization  day  and  kept  unchanged  during  the  initial¬ 
ization-assimilation  procedure.  In  this  case,  the  assimi¬ 
lated  BL  is  spread  throughout  the  model  domain  until 
the  equilibrium  is  reached  [when  the  value  of  dC/dt  is 
zero  in  Eq.  (1)].  This  provides  the  initial  BL  distribu¬ 
tion,  which  is  dynamically  balanced  with  the  physical 
conditions  at  the  time  of  the  initialization. 

According  to  (1),  the  following  equilibrium  relation 
is  reached  al  the  end  of  initialization: 
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where  C,  is  equilibrium  tracer  distribution;  subscript  i 
means  that  variables  are  taken  from  the  model  predic¬ 
tions  for  the  day  of  initialization.  The  equilibrium  field 
Cj  is  used  as  the  initial  tracer  distribution  for  the  fol¬ 


lowing  3  days’  prognostic  calculations  with  tracer  Eq. 
(1)  where  the  source-minus-sink  term  (5)  is  the  last 
term  on  the  left-hand  side  of  Eq.  (3)  [note  that  the 
source-minus-sink  term  (5)  is  not  equal  to  zero  at  the 
end  of  initialization  procedure  according  to  Eq.  (3)J. 

During  prognostic  calculations,  the  hydrodynamic 
velocities  and  diffusivities  change  in  accord  with  the 
hydrodynamic  model. 

On  the  open  boundary  of  the  model  domain,  the 
background  BL  values  (minimum  of  observed  BL)  are 
advected  into  the  model  domain  in  the  case  of  inflow, 
and  the  internal  (one  grid  inside)  BL  values  are  ad¬ 
vected  to  the  open  boundary  in  the  case  of  outflow  [see 
Shulman  et  al.  (2003)  for  more  details  about  initializa¬ 
tion  and  prognostic  simulations]. 

One  of  the  well-known  approaches  for  sensitivity 
studies  and  optimization  of  sampling  strategies  is  based 
on  the  adjoint  of  the  model  (see,  e.g..  Baker  and  Daley 
2000;  Rabier  et  al.  1996;  Errico  1997).  Also,  in  McCil- 
licuddy  et  al.  (1998),  the  adjoint  to  the  tracer  model  was 
used  for  estimation  of  the  sources  and  sinks  of  the 
population  dynamics.  Here  the  adjoint  code  to  the  ad- 
vection-diffusion-reaction  model  (1)  was  used  to  study 
the  sensitivity  of  the  modeled  BL  distributions  to  the 
sampling  strategies  of  BL  intensity. 

The  BL  forecast  measure  should  be  defined  for  a 
quantification  of  the  model  forecast  and  for  a  sensitivity 
study.  The  forecast  measure  J  can  be  any  scalar  func¬ 
tion  of  initial  conditions  such  that  its  gradient  exists.  It 
might  be  a  function  representing  a  weighted  integration 
of  C  over  a  particular  cross  section.  With  the  weight 
equaling  1 ,  the  forecast  measure  will  be  just  the  integral 
of  the  concentration.  If  the  weight  is  velocity  normal  to 
the  section,  the  forecast  measure  will  be  the  flux  of 
concentration  through  the  section.  “Intuitively,”  the 
flux  of  the  concentration  through  the  cross  section  bet¬ 
ter  represents  the  influence  of  advective  and  diffusive 
processes  on  concentration  distributions  than  just  the 
integral  of  the  concentration.  The  integrated  flux  of  the 
concentration  through  a  cross  section  in  the  bay  was 
introduced  as  forecast  measure  J  in  our  sensitivity  stud¬ 
ies. 


J  = 


Cil ,  dA  di 


(4) 


\un\  dA  dl 


where  A  is  cross  section,  it,,  is  velocity  normal  to  the 
cross  section  A,  and  T  is  time  length  of  the  forecast  ( 1-3 
days).  Note  that  J  can  be  either  positive  or  negative: 
although  C  is  positive  definite,  the  flux  can  be  of  either 
sign,  depending  on  the  orientation  of  the  velocity 
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Fig.  6.  Adjoint-based  sensitivity  maps  for  1-3  days  prior  to  the  forecast  during  the  relaxation  event  in  Aug  2000  (a)  at  section  X 
north  of  the  bay;  (b)  at  section  Y  in  the  northern  part  of  the  bay  mouth. 


through  the  section.  By  using  the  adjoint  we  estimate 
the  gradient  of  the  forecast  measure  ./  (4)  with  respect 
to  the  BL  initial  intensity  at  12-h  intervals  prior  to  the 
forecast  by 


where  s  is  the  sensitivity,  C(,  is  the  initial  three- 
dimensional  distribution  of  BL,  and  (dJ/dC{))  is  the  gra¬ 
dient  of  J  with  respect  to  initial  conditions  C0.  In  areas 
where  the  gradient  of./  has  large  positive  or  large  nega¬ 
tive  values,  a  change  in  the  BL  intensity  would  have 
created  a  large  impact  on  the  forecast.  Similarly,  in  ar¬ 
eas  where  the  gradient  is  small,  such  a  change  in  the 
initial  conditions  would  have  very  little  effect  on  the 
subsequent  forecast.  Maps  of  the  forecast  measure  gra¬ 
dients  (5)  show  the  sensitivity  of  the  forecast  to  the  BL 
distributions  prior  to  the  forecast.  Tn  other  words,  they 
show  areas  where  initial  conditions  are  most  influential 
on  the  forecast  and  where  sample  BL  intensity  is  most 
critical. 

Derivation  of  adjoint-based  sensitivity  maps  does  not 
depend  on  the  actual  BL  observations;  therefore,  opti¬ 


mal  locations  of  BL  sampling  can  be  estimated  prior  to 
the  actual  measurements. 

Here,  we  demonstrate  sensitivity  maps  for  the  relax¬ 
ation  event  during  August  2000  and  the  upwelling  event 
in  August  2003.  Maps  of  the  sensitivity  to  the  initial  BL 
distributions  for  1-3  days  prior  to  the  forecast  during 
the  relaxation  event  of  2000  are  shown  in  Figs.  6a  and 
6b  for  two  sections:  the  X  section  is  outside  of  the  bay, 
and  the  Y  section  is  inside  of  the  bay.  Large  positive 
values,  in  red,  indicate  a  significant  increase  in  the  fore¬ 
cast  measure  J  if  BL  intensity  increases  in  this  area, 
while  large  negative  values,  in  blue,  indicate  a  signifi¬ 
cant  decrease  in  the  forecast  measure  J  if  BL  intensity 
increases  in  this  area.  As  we  have  noted  before,  it  is 
important  to  sample  both  areas  where  absolute  value  of 
the  sensitivity  metric  is  high. 

Structure  in  the  positive  and  negative  regions  of  the 
sensitivity  metric  also  contains  useful  information 
about  the  propagation  of  information  within  the  sys¬ 
tem.  Current  structure  during  the  relaxation  event  (Fig. 
2)  indicates  a  predominantly  northward  flow  (at  least 
down  to  50  m)  across  section  X  (Fig.  6a).  At  the  same 
time,  the  positive  direction  of  velocity  is  defined  as 
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Flo.  7.  Adjoin t-bascd  sensitivity  maps  for  1-3  days  prior  to  the  forecast  during  the  upwelling  event  in  Aug  2003  (a)  at  section  X 
north  of  the  bay;  (b)  at  section  Y  in  the  northern  part  of  the  bay  mouth. 


southward  in  the  model.  Therefore,  any  increase  in  con¬ 
centration  C  to  the  south  of  section  X  will  contribute 
negative  values  of  product  Cu„  into  the  integrand  of 
(4),  and  therefore  decrease  the  value  of  J.  This  rela¬ 
tionship  is  why  the  gradient  [sensitivity  metric  (5)]  has 
large  negative  values  to  the  south  of  section  X.  Note 
that  the  sensitivity  maps  provide  us  with  estimates  as  to 
the  extent  of  the  area  to  which  BL  sampling  is  needed 
1 ,  2,  or  3  days  prior  to  the  forecast.  These  regions  for 
optimizing  sampling  cannot  be  deduced  using  only  the 
current  structure  alone.  For  section  Y  (Fig.  6b),  the 
direction  of  currents  has  more  variability  than  at  sec¬ 
tion  X,  and  for  this  reason  there  are  areas  of  negative 
and  positive  sensitivity.  For  section  X,  the  maps  (Fig. 
6a)  show  high  sensitivity  to  BL  intensity  in  the  northern 
part  of  the  bay  in  the  case  of  2-  and  3-day  forecasts.  This 
corresponds  to  the  BL  predictability  results  outlined  in 
Shulman  et  al.  (2003),  when  the  assimilation  of  only  the 
inside-the-bay  survey  data  into  the  tracer  model  gave  a 
good  reconstruction  of  the  observed  location  of  the  BL 
maximum  for  the  outside-the-bay  section.  At  the  same 
lime,  for  section  Y,  the  critical  areas  for  sampling  are 
inside  the  bay  and  to  the  south  of  the  bay.  Sampling  to 
the  north  and  outside  of  the  bay  is  not  important  for 


short-term  BL  predictions  inside  the  bay.  Again,  this  is 
in  agreement  with  the  results  of  short-term  BL  predict¬ 
ability  experiments  described  in  Shulman  et  al.  (2003), 
in  which  it  was  demonstrated  that  assimilation  of  ob¬ 
servations  to  the  north  of  the  bay  has  little  or  no  effect 
on  short-term  predictions  of  BL  features  observed 
along  the  section  inside  the  bay. 

Sensitivity  maps  for  the  upwelling  event  of  August 
2003  are  shown  in  Fig.  7.  For  a  section  inside  the  bay,  it 
is  critical  to  sample  inside  the  bay  and  along  the  en¬ 
trance  to  the  bay.  Sensitivity  patterns  rotate  clockwise  if 
we  move  from  L  day  to  3  days  prior  to  the  forecast.  This 
result  indicates  that  sampling  should  be  done  inside  the 
bay  counterclockwise  in  time,  which  is  supported  by  the 
presence  of  a  strong  cyclonic  eddy  inside  the  bay 
(Fig.  3). 

5.  Bioluminescence  predictability  experiments 
during  the  upwelling  event  of  August  2003 

To  test  the  sampling  strategy  indicated  by  the  sensi¬ 
tivity  map  in  Fig.  7b,  numerical  experiments  of  the  BL 
predictability  were  conducted  for  the  period  14-17  Au¬ 
gust  2003.  For  the  initialization  of  the  advection- 
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diffusion-reaction  model  (1)  on  14  August  2003,  BL 
observations  were  assimilated  into  the  model  (I),  while 
velocities  and  diffusi vities  from  frsICON  model  predic¬ 
tions  on  14  August  were  kept  unchanged  during  the 
initialization  procedure  (see  section  4).  This  procedure 
provides  the  initial  BL  distribution,  which  is  dynami¬ 
cally  balanced  with  the  physical  conditions  on  14  Au¬ 
gust.  Three  numerical  modeling  experiments  are  pre¬ 
sented  here.  Experiments  differ  in  BL  observations 
used  for  initialization  of  the  model  on  14  August.  After 
initialization,  for  all  experiments,  3  days  of  prognostic 
simulations  [when  physical  fields  in  (1)  change  in  ac¬ 
cord  with  frsICON  velocities  and  diffusivities]  were 
conducted.  Model  predictions  along  section  M  I4a  dur¬ 
ing  15-17  August  2003  were  compared  to  the  corre¬ 
sponding  BL  observations  (Fig.  5). 

a.  Experiment  1 

Only  BL  observations  from  section  Ml 4a  were  as¬ 
similated  for  initialization  of  the  model  (Fig.  4,  top). 
Therefore,  BL  intensity  only  from  the  section  of  inter¬ 
est  in  the  model  forecast  over  the  next  three  days  was 
used  in  the  initialization.  Information  from  the  adjoint 
sensitivity  map  (Fig.  7b)  was  not  used  in  this  experi¬ 
ment.  The  experiment  simulates  the  situation  that  as¬ 
sumes  that  BL  variability  is  mostly  determined  by  the 
intensity  level  of  BL  in  the  area  where  we  are  interested 
in  forecasting. 

b.  Experiment  2 

Bioluminescence  observations  from  the  three  sec¬ 
tions  M  l  4b,  H225,  and  H226  were  assimilated  for  ini¬ 
tialization  of  the  model.  This  experiment  simulates  the 
following  situation:  according  to  the  sensitivity  maps  of 
Fig.  7b,  BL  sampling  is  conducted  on  14  August  along 
three  sections,  where  the  forecast  sensitivity  is  high; 
however,  we  did  not  sample  the  area  (section  M14a) 
where  we  are  interested  in  forecasting.  In  this  case,  we 
suppose  that  BL  intensity  is  determined  mostly  by  the 
advection-diffusion  of  the  organisms  from  other  areas 
in  the  modeling  domain,  and,  therefore,  BL  intensity 
prior  to  the  forecast  sampled  in  the  area  of  interest  is 
not  very  important  for  1-3-day  forecasts. 

c.  Experiment  3 

Bioluminescence  observations  from  all  four  sections 
of  Fig.  4  (Ml4a,  M!4b,  H225.  and  H226)  were  assimi¬ 
lated  for  the  initialization  of  the  model.  This  experi¬ 
ment  simulates  the  situation  when  BL  sampling  is  con¬ 
ducted  on  14  August  along  three  sections,  where  the 
forecast  sensitivity  is  high  according  to  sensitivity  maps 


(Fig.  7b),  as  well  as  along  section  Ml 4a,  where  we  are 
interested  in  predictions.  In  this  experiment,  we  test  the 
case  when  BL  intensity  sampling  is  conducted  in  the 
area  of  the  forecast  as  well  as  in  the  most  sensitive  areas 
according  to  sensitivity  maps  of  Fig.  7b. 

Model-predicted  BL  intensity  along  section  M  14a  is 
shown  for  all  three  experiments  in  Fig.  8,  and  surface- 
predicted  model  BL  intensity  is  shown  in  Fig.  9.  These 
distributions  were  compared  to  the  observations  of  BL 
intensity  along  Ml 4a  shown  in  Fig.  5.  In  experiment  1, 
the  predicted  BL  intensity  remains  very  similar  to  the 
initial  distribution  on  14  August.  With  some  changes  in 
intensity,  BL  is  concentrated  in  the  area  around  6  km 
offshore  on  all  three  days  (Fig.  8).  The  model  in  experi¬ 
ment  l  did  not  predict  observed  offshore  translation 
and  intensification  of  BL  intensity  (Fig.  5). 

In  experiment  2  (Fig.  8),  the  model-predicted  BL  in¬ 
tensity  indicates  observed  offshore  spreading  and  inten¬ 
sification  of  BL  intensity  (Fig.  5)  along  section  M14a. 
Though  the  level  of  model-predicted  BL  intensity  is 
lower  than  the  observed  level,  the  results  of  experiment 
2  support  the  sampling  strategy  derived  from  adjoint 
sensitivity  maps  of  Fig.  7b.  Finally,  in  experiment  3  (Fig. 
8),  the  best  agreement  between  model-predicted  and 
observed  BL  intensity  is  achieved.  This  means  that  the 
B  L  sampling  according  to  the  adjoint  sensitivity  map,  as 
well  as  in  the  area  of  interest  (section  M14a),  provides 
BL  predictions  that  agree  with  the  observed  in  spatial 
and  short-term  temporal  changes. 

In  experiment  3,  Fig.  8  indicates  the  development  of 
the  secondary,  offshore  maximum  during  the  third  day 
of  predictions.  This  secondary  maximum  is  not  present 
in  the  observations  (Fig.  5).  The  second  maximum  rep¬ 
resents  the  BL  intensity  advected  from  the  north  to  the 
center  of  the  bay  (Fig.  9).  The  presence  of  this  offshore 
maximum  may  be  the  result  of  slight  deviation  of  the 
model-predicted  southward  jet  from  its  true  location. 
At  the  same  time,  in  the  case  of  longer-term  predictions 
[(9(7  days)],  the  lack  of  biological  interactions  influenc¬ 
ing  the  BL  distributions  might  well  be  a  contributor  to 
the  development  of  this  secondary  maximum. 

To  quantify  predictive  skills  in  the  experiments  con¬ 
sidered  above,  correlations  between  model-predicted 
and  observed  BL  distributions  at  section  Ml 4a  were 
calculated.  Table  I  presents  correlation  coefficients  be¬ 
tween  measurements  taken  during  14-17  August  (Fig. 
5)  and  model  initial  BL  distributions  (1.4  August,  top 
row  of  Fig.  8),  while  Table  2  presents  correlations  be¬ 
tween  observed  and  model-predicted  BL  distributions 
during  14-17  August.  Therefore,  Tables  1  and  2  provide 
comparisons  of  persistence  versus  forecast  for  each  of 
the  predictive  experiments  described  above.  Overall  for 
all  experiments,  the  forecast  has  higher  correlations 
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a.  Exp.  1  b.  Exp.  2  c.  Exp.  3 


Fi(i  8.  Model-predicted  bioluminescence  distributions  along  section  M14a.  The  distance  scale  starts  from  the  onshore  end  of 

the  section. 


with  observations  than  the  persistence.  For  example,  in 
experiment  2  (second  row  in  Table  1),  there  are  Jovv 
correlations  between  the  model  initial  and  observed  BL 
distributions  during  14-J7  August.  This  is  due  to  the 
lack  of  assimilation  of  observed  BL  at  section  M14a  on 
14  August  for  the  model  initialization  in  experiment  2. 
However,  the  forecast  in  experiment  2  (Table  2)  has 
much  higher  correlations  with  observations  than  the 
persistence  (Table  1),  especially  for  the  first  two  days. 
Note  that  the  presence  of  the  secondary  offshore  BL 
maximum  reduced  the  correlation  between  the  model 
forecast  and  observed  BL  distribution  for  experiment  3 
on  17  August. 


6.  Conclusions,  discussion,  and  future  plans 

It  has  been  demonstrated  that  the  assimilation  of  BL 
observations  into  the  adveclion-diffusion-reaction 
model  provides  a  valuable  methodology  for  the  inter¬ 
pretation  and  short-term  predictions  of  temporal  and 
spatial  changes  in  coastal  bioluminescence.  Short-term 


changes  in  observed  BL  intensity  during  wind  relax¬ 
ation  (August  2000)  and  upwelling  (August  2003) 
events  were  successfully  predicted  with  the  advection- 
diffusion-reaction  model  using  velocities  and  diffusivi- 
ties  from  the  high-resolution  circulation  model  of 
Monterey  Bay.  For  successful  short-term  BL  predic¬ 
tions,  the  optimization  of  locations  and  time  of  BL  ob¬ 
servations  prior  to  the  forecast  is  critical.  Sensitivities  of 
short-term  BL  forecasts  to  the  sampling  strategies  were 
investigated  by  using  the  adjoint  to  the  adveclion- 
diffusion-reaction  model  code.  During  the  relaxation 
event,  sensitivity  maps  of  BL  forecasts  show  that  for 
short-term  BL  intensity  predictions  at  the  cross  section 
inside  of  Monterey  Bay,  the  area  to  the  south  of  the  bay 
(Fig.  6b)  should  have  priority  when  sampling.  During 
the  upwelling  event,  sensitivity  maps  of  BL  forecasts 
suggest  a  counterclockwise  strategy  for  sampling  BL 
intensity  inside  the  bay  (Fig.  7b).  Adjoint-based  sensi¬ 
tivity  maps  provide  a  simple  methodology  for  the  opti¬ 
mization  of  BL  sampling.  Derivation  of  adjoint-based 
sensitivity  maps  does  not  depend  on  the  actual  BL  ob¬ 
servations;  therefore,  optimal  locations  of  BL  sampling 
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Fig.  9.  Surface  model-predicted  bioluminescence  distributions  in  three  numerical  experiments. 
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Taiii ,k  I.  Correlations  between  observed  and  model-predicted 
initial  BL  distributions. 


14  Aug 

15  Aug 

16  Aug 

17  Aug 

Experiment  1 

0.79 

0.48 

0.57 

0.54 

Experiment  2 

0.35 

0.33 

0.33 

0.33 

Experiment  3 

0.79 

0.63 

0.64 

0.62 

can  be  estimated  prior  to  the  actual  observations  as 
long  as  the  hydrodynamic  field  is  known  a  priori. 

Our  results  demonstrate  a  strong  dependence  of  BL 
intensity  distribution  on  flow  conditions.  During  a  re¬ 
laxation  evert,  patterns  of  cross-section  velocities  from 
the  frsICON  model  indicated  the  development  of  a 
sharp  frontal  structure  that  moved  onshore  over  3  days. 
The  observed  BL  intensity  during  these  days  also 
showed  onshore  translation,  shallowing,  and  intensifi¬ 
cation  of  the  BL  (see  Fig.  2  in  Shulman  et  al.  2003). 
Bioluminescence  observations  during  the  upwelling  of 
August  2003  (Figs.  4  and  5)  indicate  offshore  transla¬ 
tion  and  intensification  of  BL  intensity  along  section 
M  14a.  Cross-section  velocities  from  the  frsICON  model 
at  section  M  l  4a  show  offshore  translation  of  a  strong 
frontal  structure  in  the  velocity  field  (Fig.  10).  There¬ 
fore,  during  the  relaxation  event,  the  BL  intensity  maxi¬ 
mum  was  moving  onshore  with  the  velocity  front,  while 
during  the  upwelling  event  the  BL  intensity  maximum 
was  moving  offshore  with  the  velocity  front. 

The  proposed  methodology  for  short-term  BL  pre¬ 
dictions  consists  of  the  following: 

1)  setting  up  the  circulation  model  for  the  area  of  in¬ 
terest  and  conducting  24-72-h  forecasts  of  physical 
conditions; 

2)  sampling  of  BL  intensity  according  to  adjoint-based 
BL  sensitivity  maps; 

3)  dynamical  initialization  (section  4;  Shulman  et  al. 
2003)  of  the  tracer  model  when  BL  observations  are 
assimilated  into  the  tracer  model  while  velocities 
and  diffusivities  are  kept  unchanged  during  the  ini¬ 
tialization  process  (this  dynamic  initialization  proce¬ 
dure  provides  an  equilibrium  tracer  distribution  that 
is  balanced  with  the  velocity  and  diffusivity  fields 
from  the  circulation  model); 


Table  2.  Correlations  between  observed  and  model-predicted 
BL  distributions. 


14  Aug 

15  Aug 

J6  Aug 

17  Aug 

Experiment  l 

0.79 

0.63 

0.54 

0.60 

Experiment  2 

0.35 

0.54 

0.63 

0.45 

Experiment  3 

0.79 

0.67 

0.7 

0.6 

4)  using  the  equilibrium  BL  distribution  as  the  initial 
BL  field  for  3  days  of  prognostic  calculations  when 
velocities  and  diffusivities  evolve  in  accord  with  the 
circulation  model  dynamics. 

The  numerical  bioluminescence  predictability  experi¬ 
ments  demonstrated  a  close  agreement  between  ob¬ 
served  and  model-predicted  short-term  spatial  and  tem¬ 
poral  changes  of  the  coastal  bioluminescence.  Our  re¬ 
sults  also  indicate  that  advective-diffusive  processes 
largely  explain  the  short-term  evolution  of  BL  inten¬ 
sity.  At  the  same  time,  it  is  clear  that  advective- 
diffusive  processes  cannot  explain  the  full  spatial  and 
temporal  variability  of  BL  intensity,  Source  and  sink 
terms  representing  ecological  interactions,  especially 
for  relatively  long-term  predictions  [0(7  days)],  should 
be  included.  In  experiment  3,  the  secondary,  offshore, 
unobserved  maximum  developed  during  the  third  day 
of  predictions.  The  lack  of  biological  interactions  influ¬ 
encing  the  BL  distributions  might  be  at  least  a  partial 
contributor  to  the  development  of  this  secondary  maxi¬ 
mum. 

Derivation  and  parameterization  of  bioluminescence 
sources  and  sinks  represent  a  very  challenging  problem: 
the  complex  interactions  characterizing  life  cycles  of 
autotrophs,  grazers,  and  predators  producing  BL,  and, 
especially,  the  mathematical  parameterizations  and  for¬ 
mulations  of  the  biological  processes  governing  BL 
variability  in  complex  ecosystems  are  fundamental  re¬ 
search  issues.  In  McGillicuddy  et  al.  (1998)  the  sources 
and  sinks  of  the  population  dynamics  (right-hand  term 
in  the  advection-diffusion-reaction  tracer  equation) 
were  determined  by  inversion  of  the  advection-dif¬ 
fusion-reaction  tracer  equation.  Observations  from  IT 
yr  of  sampling  were  used  for  inversion,  verification,  and 
interpretation  of  population  dynamics  maps  derived 
from  the  inversion.  Our  future  research  will  include  the 
design  and  execution  of  similar  inversion  experiments 
to  determine  sources  and  sinks  of  BL  intensity  in 
Monterey  Bay.  The  adjoint  model  developed  for  the 
ICON  model  tracer  routine  will  be  the  basis  for  these 
inversion  experiments.  The  results  of  these  inversion 
experiments,  together  with  research  presented  here, 
will  provide  an  understanding  of  the  relative  contribu¬ 
tions  of  advective-diffusive  versus  biological  processes 
to  short-term  and  relatively  long-term  [0(7  days)]  vari¬ 
ability  of  BL  intensity  in  the  bay.  An  important  tool  for 
achieving  this  level  of  understanding  will  be  the 
Coupled  Physical  BioOptical  Model  (http://www7320. 
nrlssc.navy.mil/cobiopp/),  as  used  for  the  interpretation 
of  bioluminescence  inversion  experiments  and  for  the 
development  of  numerical  parameterizations  of  biolu¬ 
minescence  sources  and  sinks. 
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14th  CroM-8*ctlor»l  Vtlodty  (m/i)  |  Y  |  Day  +0 


Distarco  from  Slartol  Sccion  (lun) 


DatBrc«  from  Start ot  Sccton(km) 


Distance  from  Start  ot  Sectt>n  (km) 


Detarcefrom  Slarto)  Section  (kin) 


Fid.  10.  ICON-model-prcdicied  cross-seciion  velocities  along  section  Y.  Positive  cross- 
sectional  velocity  is  from  southwest  to  northeast  across  section  Y.  The  distance  scale  starts 
from  the  onshore  end  of  the  section. 
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